package com.qyx.pmpucat.util;

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.integration.TrapezoidIntegrator;
import org.apache.commons.math3.distribution.NormalDistribution;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.RealMatrix;

public class EapSem {

    private static final double D = 1.0; // 这个值可以在构造函数中传入

    public static class PiResult {
        public double[][] Pi;  // 概率矩阵
        public double[][] dPi;  // 一阶导数矩阵
        public double[][] d2Pi;  // 二阶导数矩阵
        public double[][] d3Pi;  // 三阶导数矩阵

        // 构造函数，初始化各个矩阵
        public PiResult(double[][] Pi, double[][] dPi, double[][] d2Pi, double[][] d3Pi) {
            this.Pi = Pi;
            this.dPi = dPi;
            this.d2Pi = d2Pi;
            this.d3Pi = d3Pi;
        }
    }

    /**
     * 计算项目反应理论 (IRT) 中的 Pi 矩阵及其导数
     * 该方法主要用于计算在给定能力值（theta）、项目参数和缩放因子D的情况下，每个项目的概率矩阵Pi及其一阶、二阶和三阶导数
     *
     * @param th 被试的能力参数
     * @param it 项目参数矩阵，每行对应一个项目，第一列是区分度参数a，后续列是难度参数b
     * @param D 缩放因子，通常为1.7，用于调整能力尺度
     * @return 返回一个PiResult对象，包含Pi矩阵及其导数
     */
    public static PiResult Pi(double th, double[][] it, double D) {
        // 将项目参数数组转换为RealMatrix对象，便于后续操作
        RealMatrix itemParams = new Array2DRowRealMatrix(it);
        // 获取项目的数量和最大类别数
        int nItems = itemParams.getRowDimension();
        int maxCategories = itemParams.getColumnDimension();

        // 初始化概率矩阵及其一阶、二阶和三阶导数矩阵
        double[][] prov = new double[nItems][maxCategories + 1];
        double[][] prov1 = new double[nItems][maxCategories + 1];
        double[][] prov2 = new double[nItems][maxCategories + 1];
        double[][] prov3 = new double[nItems][maxCategories + 1];

        // 遍历每个项目
        for (int i = 0; i < nItems; i++) {
            // 提取当前项目的区分度参数a和难度参数b数组
            double aj = itemParams.getEntry(i, 0);
            double[] bj = itemParams.getSubMatrix(i, i, 1, maxCategories - 1).getData()[0];
            // 初始化概率Pj及其导数数组
            double[] Pjs = new double[bj.length + 2];
            Pjs[0] = 1.0;
            // 计算每个类别概率Pj
            for (int j = 0; j < bj.length; j++) {
                double ej = Math.exp(D * aj * (th - bj[j]));
                Pjs[j + 1] = ej / (1 + ej);
            }
            Pjs[Pjs.length - 1] = 0.0;

            // 初始化并计算概率Pj的一阶、二阶和三阶导数
            double[] dPjs = new double[Pjs.length];
            double[] d2Pjs = new double[Pjs.length];
            double[] d3Pjs = new double[Pjs.length];

            for (int j = 0; j < Pjs.length; j++) {
                dPjs[j] = D * aj * Pjs[j] * (1 - Pjs[j]);
                d2Pjs[j] = D * aj * (dPjs[j] - 2 * Pjs[j] * dPjs[j]);
                d3Pjs[j] = D * aj * (d2Pjs[j] - 2 * Math.pow(dPjs[j], 2) - 2 * Pjs[j] * d2Pjs[j]);
            }

            // 计算并填充概率矩阵Pi及其导数矩阵
            for (int j = 0; j < maxCategories; j++) {
                prov[i][j] = Pjs[j] - Pjs[j + 1];
                prov1[i][j] = dPjs[j] - dPjs[j + 1];
                prov2[i][j] = d2Pjs[j] - d2Pjs[j + 1];
                prov3[i][j] = d3Pjs[j] - d3Pjs[j + 1];
            }
        }

        // 返回计算结果
        return new PiResult(prov, prov1, prov2, prov3);
    }

    /**
     * 使用梯形积分法进行数值积分
     *
     * @param x X轴上的数据点数组，用于插值
     * @param y 对应x数据点的Y轴值数组
     * @return 返回积分的结果
     *
     * 此方法采用梯形积分法来逼近由离散数据点定义的曲线下的面积
     * 它首先创建一个梯形积分器实例，然后定义一个基于给定数据点的插值函数
     * 最后，使用梯形积分法在给定的区间上对这个插值函数进行积分
     */
    private static double integrateCat(double[] x, double[] y) {
        // 创建梯形积分器实例
        TrapezoidIntegrator integrator = new TrapezoidIntegrator();
        // 定义插值函数，基于给定的数据点
        UnivariateFunction function = t -> interpolate(x, y, t);
        // 执行积分操作，参数分别为：积分的最大次数，被积函数，积分区间的下限和上限
        return integrator.integrate(2000, function, x[0], x[x.length - 1]);
    }

    /**
     * 线性插值函数
     * 该函数用于根据给定的x和y数组以及一个目标值t，计算出对应的插值y值
     * 线性插值是通过找到目标值t在x数组中的两个相邻点，然后根据这两个点的y值进行线性插值计算
     *
     * @param x 一个有序的double数组，代表x轴上的已知数据点
     * @param y 一个与x数组相对应的double数组，代表y轴上的已知数据点
     * @param t 目标值，即需要进行插值的x轴上的点
     * @return 返回插值计算得到的y轴上的值
     * @throws IllegalArgumentException 如果目标值t超出了x数组的范围，则抛出此异常
     */
    private static double interpolate(double[] x, double[] y, double t) {
        // 如果目标值t小于等于x数组的第一个元素，返回y数组的第一个元素
        if (t <= x[0]) return y[0];
        // 如果目标值t大于等于x数组的最后一个元素，返回y数组的最后一个元素
        if (t >= x[x.length - 1]) return y[y.length - 1];
        // 遍历x数组，寻找目标值t所在两个相邻点之间的位置
        for (int i = 0; i < x.length - 1; i++) {
            // 当找到目标值t所在的位置区间时，进行线性插值计算
            if (x[i] <= t && t <= x[i + 1]) {
                // 线性插值公式：y = y[i] + (y[i + 1] - y[i]) * (t - x[i]) / (x[i + 1] - x[i])
                return y[i] + (y[i + 1] - y[i]) * (t - x[i]) / (x[i + 1] - x[i]);
            }
        }
        // 如果目标值t超出x数组的范围，抛出异常
        throw new IllegalArgumentException("Interpolation point out of range.");
    }

    /**
     * 计算给定参数下的预期后验标准误（EAP SEM）
     * 该方法用于估计受试者能力参数的精确度
     *
     * @param thEst 受试者的估计能力参数
     * @param it 项目特征曲线参数的二维数组，每行对应一个项目，每列依次表示项目区分度、难度参数
     * @param x 受试者在每个项目上的反应模式
     * @return 返回受试者能力参数的预期后验标准误
     */
    public static double eapSem1(double thEst, double[][] it, int[] x) {
        // 定义对数似然函数 LL
        UnivariateFunction LL = th -> {
            // 计算每个项目在当前能力参数下的概率
            PiResult piRes = Pi(th, it, D);
            double res = 1.0;
            // 计算每个项目的似然值
            for (int i = 0; i < x.length; i++) {
                res *= piRes.Pi[i][x[i]];
            }
            return res;
        };
        // 创建梯形积分器实例
        TrapezoidIntegrator integrator = new TrapezoidIntegrator();
        // 生成积分用的点序列
        double[] X = linspace(-4, 4, 33);

        // 计算分子部分的积分
        double numeratorIntegral = integrate(s -> gg(s, thEst, LL), X, integrator);
        // 计算分母部分的积分
        double denominatorIntegral = integrate(s -> hh(s, LL), X, integrator);

        // 返回预期后验标准误
        return Math.sqrt(numeratorIntegral / denominatorIntegral);
    }

    /**
     * 计算给定参数下的似然值
     * 该方法用于计算在特定参数条件下，观测数据出现的概率
     * 它通过调用Pi方法计算每个观测值的概率，并将这些概率相乘得到最终的似然值
     *
     * @param th 代表特定的参数值
     * @param it 包含多个观测值的数组，每个观测值可能有多个属性
     * @param x 观测结果的数组，每个元素代表对应观测值的结果
     * @return 返回计算得到的似然值
     */
    private static double likelihood(double th, double[][] it, int[] x) {
        // 调用Pi方法，计算在参数th下，每个观测值的概率
        PiResult piResult = Pi(th, it, D);
        // 初始化概率值为1，准备进行乘法运算
        double prob = 1.0;
        // 遍历观测结果数组
        for (int i = 0; i < x.length; i++) {
            // 将每个观测值的概率乘以当前概率值，得到最终的似然值
            prob *= piResult.Pi[i][x[i] + 1];
        }
        // 返回计算得到的似然值
        return prob;
    }

    private static double gg(double s, double thEst, UnivariateFunction LL) {
        NormalDistribution normalDist = new NormalDistribution(0, 1);
        return Math.pow(s - thEst, 2) * normalDist.density(s) * LL.value(s);
    }

    private static double hh(double s, UnivariateFunction LL) {
        NormalDistribution normalDist = new NormalDistribution(0, 1);
        return normalDist.density(s) * LL.value(s);
    }

    /**
     * 生成一个等差数列
     * 该方法用于在指定的起始值和结束值之间，生成一个包含指定数量元素的等差数列
     *
     * @param start 起始值，数列的第一个元素
     * @param end 结束值，数列的最后一个元素
     * @param num 数列中元素的数量
     * @return 返回一个包含等差数列的数组
     */
    private static double[] linspace(double start, double end, int num) {
        // 创建一个数组来存储等差数列
        double[] sequence = new double[num];
        // 计算等差数列的公差
        double step = (end - start) / (num - 1);
        // 使用循环来填充数组，生成等差数列
        for (int i = 0; i < num; i++) {
            // 计算数列中的每个元素，并将其存储到数组中
            sequence[i] = start + i * step;
        }
        // 返回生成的等差数列数组
        return sequence;
    }

    /**
     * 使用梯形法则近似计算给定函数在一系列点范围内的积分
     *
     * @param f 要积分的单变量函数
     * @param points 定义积分区间的点数组，积分将在这些点定义的子区间上进行
     * @param integrator 实现了梯形积分法的对象，用于执行实际的积分计算
     * @return 在所有子区间上计算得到的积分总和
     */
    private static double integrate(UnivariateFunction f, double[] points, TrapezoidIntegrator integrator) {
        double integral = 0.0;
        // 遍历所有子区间，计算每个子区间的积分并累加
        for (int i = 0; i < points.length - 1; i++) {
            // 调用integrator的integrate方法计算每个子区间的积分，结果累加到integral变量中
            integral += integrator.integrate(10000, f, points[i], points[i + 1]);
        }
        // 返回累加的积分值作为最终结果
        return integral;
    }


}